library(lme4)
## Loading required package: Matrix
library(nlme)
##
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
##
## lmList
library(car)
## Loading required package: carData
library(ggplot2)
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:nlme':
##
## collapse
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(performance)
library(emmeans)
library(devtools)
## Loading required package: usethis
##
## Attaching package: 'devtools'
## The following object is masked from 'package:emmeans':
##
## test
devtools::install_github("psyteachr/introdataviz")
## Using GitHub PAT from the git credential store.
## Skipping install of 'introdataviz' from a github remote, the SHA1 (0519c980) has not changed since last install.
## Use `force = TRUE` to force installation
datastore <- read.csv("~/Temnos2022Venom/dissectiondataset.csv")
data<-na.omit(datastore)
summarydata<-read.csv("~/Temnos2022Venom/summaryraiddata2022.csv")
Data set has both workers and queens. Here are the variables in the dataset:
Date.Dissected: Date ant was dissected MM/DD/YY.
Date.Frozen: Date ant was frozen for dissection. This is approximately 2 weeks after raid trial ended MM/DD/YY.
Date.Frozen: Date colony was collected MM/DD/YY.
Collection.Site: Site at which nest containing ant was collected.
Nest: Number identifier for nest.
Pair.Group: Number identifier for pair group. Each colony collected in the field was divided between 2 nests as a matched pair. Nests with the same Pair group are genetically related (came from the same colony) but one nest in each pair experienced a raiding trial, and one did not.
Raided: Whether or not the nest was raided. Y means yes, N means no.
Worker.Queen : Whether the dissected ant was a worker or a queen
Behavior.Code: If Worker.Queen==” Worker”, this is N for nurse, F for forager, or X for unclassified. If Worker.Queen == “Queen”, this is AQ for alate queen or DQ for dealate queen.
Ant : Ant number (identifier)
Webers.Length.mm : Weber’s length in mm
Venom.Sac.Length.mm : Length of venom sac in mm
Venom.Sac.Width.mm : Width of venom sac in mm
Final.Worker.Count : The number of workers in the nest on the day it was censused/frozen.
Final.Queen.Count.Dealate : The number of queens without wings in the nest on the day it was censused/frozen.
Final.Pupae.Count : The number of pupae in the nest on the day it was censused/frozen.
Final.Larvae.Count : The number of larvae in the nest on the day it was censused/frozen.
Final.Egg.Count : The approximate number of eggs in the nest on the day it was censused/frozen.
Final.Queen.Count.Alate : The number of winged queens in the nest on the day it was censused/frozen.
Final.Male.Count: The number of males in the nest on the day it was censused/frozen.
In this chunk, I add a few new variables to the dataset that can be calculated from the initial variables. I also make sure the date is in the correct format. The most important variable created here is Venom.Volume, which is calculated with the formula Venom.Volume \(=\frac{\pi}{6} \times(L \times W^2)\) for each venom sac. I also create a dataset that is only workers (subset of dataset)
data$Raided<-as.factor(data$Raided)
data$Date.Dissected <- as.Date((data$Date.Dissected), format = "%m/%d/%y")
data$Date.Frozen <- as.Date((data$Date.Frozen), format = "%m/%d/%y")
data$Date.Collected <- as.Date((data$Date.Collected), format = "%m/%d/%y")
#In R, first day is automatically Jan 1, 1970. This means that when you put a fixed effect in date format in a model, the results become super difficule to interpret. Therefore we must create a new date variable to include in the model to make results more interpretable. To do this I first convert Frozen date to a numeric, where day 1 is the first day I froze ant nests, August 24, 2022 (since first day in R is 1970 I convert to a numeric and then subtract every day before that, which is 19226 days)
data$Day.Frozen<-as.numeric(data$Date.Frozen)-19226
data$Final.Brood=data$Final.Larvae.Count+data$Final.Pupae.Count
data$Final.Brood.Worker.Ratio=data$Final.Brood/data$Final.Worker.Count
#Scale variables to make results more interpretable based on summary data. Meaning because every ant within a colony has the same final worker count and final pupae count, we need to scale based on the colony level and then give those numbers to the individual level dataset
meanworkers<-mean(summarydata$Final.Worker.Count)
SDworkers<-sd(summarydata$Final.Worker.Count)
data$Worker.Count.Scaled<-(data$Final.Worker.Count-meanworkers)/SDworkers
meanpupae<-mean(summarydata$Final.Pupae.Count)
SDpupae<-sd(summarydata$Final.Pupae.Count)
data$Pupae.Count.Scaled<-(data$Final.Pupae.Count-meanpupae)/SDpupae
meanlarvae<-mean(summarydata$Final.Larvae.Count)
SDlarvae<-sd(summarydata$Final.Larvae.Count)
data$Larvae.Count.Scaled<-(data$Final.Larvae.Count-meanlarvae)/SDlarvae
data$Collection.Site<-as.factor(data$Collection.Site)
#Create a variable for venom volume, which is an ellipsoid using venom sac length and width. Since both length and width are measured in mm, this volume is in mm^3, which is 1:1 with microliters. For the purposes of graphing and modeling this problem, we will multiply by 1000 to convert to nanoliters
data$Venom.Volume<-(pi/6)*data$Venom.Sac.Length.mm*data$Venom.Sac.Width.mm^2*1000
#Create a dataset for just workers
workerdata<-subset(data,Worker.Queen=="Worker")
queendata<-subset(data,Worker.Queen=="Queen")
To investigate this we test whether behavior code can predict weber’s length in our data set. We ran two models, one with all ants in the data set and one with just ants from unraided colonies in case ants died during a raid and more ants in the colony switched castes. Regardless, behavioral caste is not a statistically significant predictor of Weber’s length.
sizemodel<-lme(Webers.Length.mm~Behavior.Code,data=workerdata, random=~1|Nest)
summary(sizemodel)
## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## -2062.388 -2040.581 1036.194
##
## Random effects:
## Formula: ~1 | Nest
## (Intercept) Residual
## StdDev: 0.03227909 0.03746912
##
## Fixed effects: Webers.Length.mm ~ Behavior.Code
## Value Std.Error DF t-value p-value
## (Intercept) 0.7968349 0.007648137 552 104.18680 0.0000
## Behavior.CodeN 0.0077771 0.005444028 552 1.42855 0.1537
## Behavior.CodeX 0.0056124 0.005079065 552 1.10500 0.2696
## Correlation:
## (Intr) Bhv.CN
## Behavior.CodeN -0.509
## Behavior.CodeX -0.537 0.764
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.59849431 -0.59716620 0.03381316 0.64082554 3.01559231
##
## Number of Observations: 582
## Number of Groups: 28
plot(sizemodel)
qqnorm(resid(sizemodel))
hist(resid(sizemodel))
sizemodelunraidedonly<-lme(Webers.Length.mm~Behavior.Code,data=subset(workerdata,Raided=="N"), random=~1|Nest)
summary(sizemodelunraidedonly)
## Linear mixed-effects model fit by REML
## Data: subset(workerdata, Raided == "N")
## AIC BIC logLik
## -1011.311 -992.9269 510.6553
##
## Random effects:
## Formula: ~1 | Nest
## (Intercept) Residual
## StdDev: 0.03416408 0.0387244
##
## Fixed effects: Webers.Length.mm ~ Behavior.Code
## Value Std.Error DF t-value p-value
## (Intercept) 0.7934376 0.011310587 279 70.15000 0.0000
## Behavior.CodeN 0.0055531 0.007816701 279 0.71041 0.4780
## Behavior.CodeX 0.0105755 0.007420036 279 1.42526 0.1552
## Correlation:
## (Intr) Bhv.CN
## Behavior.CodeN -0.503
## Behavior.CodeX -0.525 0.775
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.32184914 -0.59486995 0.05129671 0.62540029 3.06965088
##
## Number of Observations: 295
## Number of Groups: 14
plot(sizemodelunraidedonly)
qqnorm(resid(sizemodelunraidedonly))
hist(resid(sizemodelunraidedonly))
violinbynest<-ggplot(data, aes(x=as.factor(Nest), y=Venom.Volume))+geom_violin(aes(color=Raided)) + geom_boxplot(width=0.1)
violinbynest
violinraided<-ggplot(subset(data,Worker.Queen=="Worker"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraided
violinraidedforager<-ggplot(subset(workerdata,Behavior.Code=="F"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraidedforager
violinraidednurse<-ggplot(subset(workerdata,Behavior.Code=="N"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraidednurse
violinraidedother<-ggplot(subset(workerdata,Behavior.Code=="X"), aes(x=Raided, y=Venom.Volume))+geom_violin()+ geom_boxplot(width=0.1)
violinraidedother
#The way Raid.Payoff is currently in dataset is that the raid payoff for the raided colony is also in the corresponding unraided matched pair. This does not make sense to include in the model, so this loop goes through each line and says if the colony didn't experience a raid, the payoff was zero.
for (i in 1:(nrow(workerdata))){
if (workerdata$Raided[i]=="N"){workerdata$Raid.Payoff.Total[i]=0}
}
workerdata$Final.Pupae.Worker.Ratio<-workerdata$Final.Pupae.Count/workerdata$Final.Worker.Count
workerdata$Final.Brood.And.Eggs<-workerdata$Final.Brood+workerdata$Final.Egg.Count
model<-lme(Venom.Volume~Raided*Behavior.Code+Day.Frozen+ Final.Brood.And.Eggs+Final.Worker.Count+Final.Pupae.Worker.Ratio+I(Webers.Length.mm^3), data=workerdata, random=list(Collection.Site=~1,Pair.Group=~1,Nest=~1))
summary(model)
## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 718.4553 783.6661 -344.2276
##
## Random effects:
## Formula: ~1 | Collection.Site
## (Intercept)
## StdDev: 0.02427888
##
## Formula: ~1 | Pair.Group %in% Collection.Site
## (Intercept)
## StdDev: 0.04567579
##
## Formula: ~1 | Nest %in% Pair.Group %in% Collection.Site
## (Intercept) Residual
## StdDev: 0.07610764 0.4142232
##
## Fixed effects: Venom.Volume ~ Raided * Behavior.Code + Day.Frozen + Final.Brood.And.Eggs + Final.Worker.Count + Final.Pupae.Worker.Ratio + I(Webers.Length.mm^3)
## Value Std.Error DF t-value p-value
## (Intercept) 0.1768862 0.1879601 549 0.941084 0.3471
## RaidedY 0.2914069 0.1051427 10 2.771538 0.0197
## Behavior.CodeN -0.1870605 0.0828191 549 -2.258664 0.0243
## Behavior.CodeX -0.1359979 0.0788279 549 -1.725251 0.0850
## Day.Frozen 0.0087471 0.0043938 8 1.990810 0.0817
## Final.Brood.And.Eggs -0.0008489 0.0011082 10 -0.766080 0.4613
## Final.Worker.Count 0.0017806 0.0025540 10 0.697175 0.5016
## Final.Pupae.Worker.Ratio -1.4126389 1.0496691 10 -1.345794 0.2081
## I(Webers.Length.mm^3) 1.6375319 0.2242739 549 7.301483 0.0000
## RaidedY:Behavior.CodeN -0.2640293 0.1193151 549 -2.212874 0.0273
## RaidedY:Behavior.CodeX -0.2406044 0.1111244 549 -2.165180 0.0308
## Correlation:
## (Intr) RaiddY Bhv.CN Bhv.CX Dy.Frz F.B.A. Fn.W.C
## RaidedY -0.282
## Behavior.CodeN -0.279 0.575
## Behavior.CodeX -0.262 0.604 0.775
## Day.Frozen -0.532 0.025 -0.014 -0.022
## Final.Brood.And.Eggs -0.328 0.038 -0.033 -0.007 0.731
## Final.Worker.Count -0.445 0.014 -0.012 -0.066 -0.045 -0.259
## Final.Pupae.Worker.Ratio -0.281 0.069 0.051 0.057 0.184 0.134 0.237
## I(Webers.Length.mm^3) -0.640 -0.032 -0.044 -0.071 -0.003 -0.165 0.228
## RaidedY:Behavior.CodeN 0.199 -0.811 -0.692 -0.538 0.050 0.029 0.015
## RaidedY:Behavior.CodeX 0.201 -0.866 -0.548 -0.706 0.009 -0.008 0.026
## F.P.W. I(W.L. RY:B.CN
## RaidedY
## Behavior.CodeN
## Behavior.CodeX
## Day.Frozen
## Final.Brood.And.Eggs
## Final.Worker.Count
## Final.Pupae.Worker.Ratio
## I(Webers.Length.mm^3) 0.040
## RaidedY:Behavior.CodeN -0.049 -0.008
## RaidedY:Behavior.CodeX -0.040 0.047 0.767
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.06580499 -0.74442354 -0.06137433 0.53783842 3.76752120
##
## Number of Observations: 582
## Number of Groups:
## Collection.Site
## 5
## Pair.Group %in% Collection.Site
## 14
## Nest %in% Pair.Group %in% Collection.Site
## 28
emmeans(model, revpairwise ~ Raided|Behavior.Code)
## $emmeans
## Behavior.Code = F:
## Raided emmean SE df lower.CL upper.CL
## N 1.141 0.0768 4 0.928 1.35
## Y 1.433 0.0781 4 1.216 1.65
##
## Behavior.Code = N:
## Raided emmean SE df lower.CL upper.CL
## N 0.954 0.0522 4 0.809 1.10
## Y 0.982 0.0563 4 0.825 1.14
##
## Behavior.Code = X:
## Raided emmean SE df lower.CL upper.CL
## N 1.005 0.0448 4 0.881 1.13
## Y 1.056 0.0447 4 0.932 1.18
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## Behavior.Code = F:
## contrast estimate SE df t.ratio p.value
## Y - N 0.2914 0.1051 10 2.772 0.0197
##
## Behavior.Code = N:
## contrast estimate SE df t.ratio p.value
## Y - N 0.0274 0.0703 10 0.389 0.7051
##
## Behavior.Code = X:
## contrast estimate SE df t.ratio p.value
## Y - N 0.0508 0.0562 10 0.904 0.3871
##
## Degrees-of-freedom method: containment
joint_tests(model)
## model term df1 df2 F.ratio p.value
## Raided 1 10 5.545 0.0403
## Behavior.Code 2 549 14.541 <.0001
## Day.Frozen 1 8 3.963 0.0817
## Final.Brood.And.Eggs 1 10 0.587 0.4613
## Final.Worker.Count 1 10 0.486 0.5016
## Final.Pupae.Worker.Ratio 1 10 1.811 0.2081
## Webers.Length.mm 1 549 53.312 <.0001
## Raided:Behavior.Code 2 549 2.714 0.0672
plot(model)
qqnorm(resid(model))
hist(resid(model))
#Controlling for body size What is the correct way to control for body size? This chunk runs the same model run in the previous chunk but controls for weber’s length cubed (a measurement of body size) by dividing the explanatory variable by it as opposed to using it as a fixed effect. The results are qualitatively the same and the fit is not as good as in the previous model.
model2<-lme(Venom.Volume/I(Webers.Length.mm^3)~Raided*Behavior.Code+Day.Frozen+ Final.Brood.And.Eggs+Final.Worker.Count+Final.Pupae.Worker.Ratio, data=workerdata, random=list(Collection.Site=~1,Pair.Group=~1,Nest=~1))
summary(model2)
## Linear mixed-effects model fit by REML
## Data: workerdata
## AIC BIC logLik
## 1498.201 1559.089 -735.1004
##
## Random effects:
## Formula: ~1 | Collection.Site
## (Intercept)
## StdDev: 0.05681686
##
## Formula: ~1 | Pair.Group %in% Collection.Site
## (Intercept)
## StdDev: 0.0001041914
##
## Formula: ~1 | Nest %in% Pair.Group %in% Collection.Site
## (Intercept) Residual
## StdDev: 0.1529016 0.8219203
##
## Fixed effects: Venom.Volume/I(Webers.Length.mm^3) ~ Raided * Behavior.Code + Day.Frozen + Final.Brood.And.Eggs + Final.Worker.Count + Final.Pupae.Worker.Ratio
## Value Std.Error DF t-value p-value
## (Intercept) 2.040822 0.2745205 550 7.434132 0.0000
## RaidedY 0.526210 0.2085407 10 2.523297 0.0302
## Behavior.CodeN -0.350172 0.1639718 550 -2.135566 0.0332
## Behavior.CodeX -0.268830 0.1558469 550 -1.724966 0.0851
## Day.Frozen 0.014603 0.0082171 8 1.777152 0.1134
## Final.Brood.And.Eggs -0.002750 0.0020377 10 -1.349605 0.2069
## Final.Worker.Count 0.005464 0.0044681 10 1.222963 0.2494
## Final.Pupae.Worker.Ratio -3.658864 2.0126765 10 -1.817909 0.0991
## RaidedY:Behavior.CodeN -0.509998 0.2366036 550 -2.155497 0.0316
## RaidedY:Behavior.CodeX -0.428976 0.2200864 550 -1.949125 0.0518
## Correlation:
## (Intr) RaiddY Bhv.CN Bhv.CX Dy.Frz F.B.A. Fn.W.C
## RaidedY -0.403
## Behavior.CodeN -0.416 0.574
## Behavior.CodeX -0.418 0.604 0.775
## Day.Frozen -0.693 0.021 -0.017 -0.024
## Final.Brood.And.Eggs -0.586 0.029 -0.043 -0.020 0.760
## Final.Worker.Count -0.373 0.017 0.000 -0.055 -0.055 -0.226
## Final.Pupae.Worker.Ratio -0.331 0.066 0.052 0.059 0.181 0.134 0.258
## RaidedY:Behavior.CodeN 0.260 -0.811 -0.693 -0.540 0.056 0.033 0.018
## RaidedY:Behavior.CodeX 0.312 -0.866 -0.547 -0.705 0.012 0.002 0.016
## F.P.W. RY:B.CN
## RaidedY
## Behavior.CodeN
## Behavior.CodeX
## Day.Frozen
## Final.Brood.And.Eggs
## Final.Worker.Count
## Final.Pupae.Worker.Ratio
## RaidedY:Behavior.CodeN -0.046
## RaidedY:Behavior.CodeX -0.040 0.768
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.58930630 -0.70688091 -0.07944495 0.55085331 5.06421913
##
## Number of Observations: 582
## Number of Groups:
## Collection.Site
## 5
## Pair.Group %in% Collection.Site
## 14
## Nest %in% Pair.Group %in% Collection.Site
## 28
emmeans(model2, revpairwise ~ Raided|Behavior.Code)
## $emmeans
## Behavior.Code = F:
## Raided emmean SE df lower.CL upper.CL
## N 2.28 0.1510 4 1.86 2.70
## Y 2.81 0.1537 4 2.38 3.23
##
## Behavior.Code = N:
## Raided emmean SE df lower.CL upper.CL
## N 1.93 0.1024 4 1.65 2.21
## Y 1.95 0.1096 4 1.64 2.25
##
## Behavior.Code = X:
## Raided emmean SE df lower.CL upper.CL
## N 2.01 0.0872 4 1.77 2.25
## Y 2.11 0.0863 4 1.87 2.35
##
## Degrees-of-freedom method: containment
## Results are given on the identity (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## Behavior.Code = F:
## contrast estimate SE df t.ratio p.value
## Y - N 0.5262 0.209 10 2.523 0.0302
##
## Behavior.Code = N:
## contrast estimate SE df t.ratio p.value
## Y - N 0.0162 0.139 10 0.116 0.9098
##
## Behavior.Code = X:
## contrast estimate SE df t.ratio p.value
## Y - N 0.0972 0.112 10 0.871 0.4041
##
## Note: contrasts are still on the identity scale
## Degrees-of-freedom method: containment
joint_tests(model2)
## model term df1 df2 F.ratio p.value
## Raided 1 10 4.203 0.0675
## Behavior.Code 2 550 13.354 <.0001
## Day.Frozen 1 8 3.158 0.1134
## Final.Brood.And.Eggs 1 10 1.821 0.2069
## Final.Worker.Count 1 10 1.496 0.2494
## Final.Pupae.Worker.Ratio 1 10 3.305 0.0991
## Raided:Behavior.Code 2 550 2.428 0.0892
plot(model2)
qqnorm(resid(model2))
hist(resid(model2))
Do nurses and foragers have the same amount of venom? Graph it.
venombybehavior<-ggplot(workerdata,aes(Behavior.Code,Venom.Volume,fill=Raided))+
geom_boxplot()+geom_point(alpha=0.4,position = position_jitterdodge(jitter.width=0.3))+
xlab("Behavioral Group")+ylab("Venom Volume (nL)")+theme_bw(base_size = 22)+
scale_x_discrete(labels=c('Foragers', 'Nurses', 'Unclassified'))+scale_fill_brewer(palette = "Set2",labels = c("No", "Yes"))
venombybehavior
ggsave('venombybehavior.svg',plot=venombybehavior,dpi=900,width =35, height = 20, units = "cm",device='svg')
venombybehavior<-ggplot(workerdata,aes(Behavior.Code,Venom.Volume,fill=Raided))+
introdataviz::geom_split_violin(alpha = .4, trim = FALSE) +
geom_boxplot(width = .2, alpha = .6, fatten = NULL, show.legend = FALSE) +
stat_summary(fun.data = "mean_se", geom = "pointrange", show.legend = F,
position = position_dodge(.175)) +
xlab("Behavioral Group")+scale_y_continuous(name = "Venom Volume (nL)",
breaks = seq(0, 3, 0.5),
limits = c(0, 3.5)) +theme_bw(base_size = 18)+
scale_x_discrete(labels=c('Foragers', 'Nurses', 'Unclassified'))+scale_fill_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes")) +
theme_minimal()
venombybehavior
## Warning: Removed 233 rows containing missing values or values outside the scale range
## (`geom_split_violin()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
ggsave('venombybehavior.svg',plot=venombybehavior,dpi=900,width =20, height = 12, units = "cm",device='svg')
## Warning: Removed 233 rows containing missing values or values outside the scale range
## (`geom_split_violin()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
#Getting warning about NAs, but there are not any
sum(is.na(workerdata$Behavior.Code))
## [1] 0
sum(is.na(workerdata$Venom.Volume))
## [1] 0
sum(is.na(workerdata$Raided))
## [1] 0
#COUNT THE ANTS
raided<-subset(workerdata,Raided=="Y")
unraided<-subset(workerdata,Raided=="N")
nrow(subset(workerdata,Behavior.Code=="N")) #179 nurses
## [1] 179
nrow(subset(raided,Behavior.Code=="N")) #82 raided, 97 unraided
## [1] 82
nrow(subset(unraided,Behavior.Code=="N"))
## [1] 97
nrow(subset(workerdata,Behavior.Code=="F")) #69 foragers
## [1] 69
nrow(subset(raided,Behavior.Code=="F")) #34 raided, 35 unraided
## [1] 34
nrow(subset(unraided,Behavior.Code=="F"))
## [1] 35
nrow(subset(workerdata,Behavior.Code=="X")) #334 unclassified
## [1] 334
nrow(subset(raided,Behavior.Code=="X")) #171 raided, 163 unraided
## [1] 171
nrow(subset(unraided,Behavior.Code=="X"))
## [1] 163
#foragers
ggplot(subset(workerdata, Behavior.Code=="F"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +geom_point()+geom_smooth(method = "lm", fill = NA)
## `geom_smooth()` using formula = 'y ~ x'
ggplot(subset(workerdata, Behavior.Code=="N"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +geom_point()+geom_smooth(method = "lm", fill = NA)
## `geom_smooth()` using formula = 'y ~ x'
Graphical depiction of venom volume by weber’s length
#Graphs by Weber's Length
#overall
weblengthworkers<-ggplot(workerdata, aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided)) +
geom_point()+geom_smooth(method = "lm", fill = NA)+theme_bw(base_size = 22)+scale_y_continuous(name="Venom volume (nL)", breaks = seq(0, 3, 0.5))+scale_x_continuous(name = ("Weber's length (mm)"),breaks = seq(0.6,1,0.1))+scale_color_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes"))+ggtitle("A")+theme(plot.margin = margin(1, 1, 1, 1, "cm"))
ggsave('weblengthworkers.svg',plot=weblengthworkers,dpi=900,width =25, height = 20, units = "cm",device='svg')
## `geom_smooth()` using formula = 'y ~ x'
weblengthforagers<-ggplot(subset(workerdata,Behavior.Code=="F"), aes(x=Webers.Length.mm, y=Venom.Volume, group=Raided, color=Raided))+ geom_point()+geom_smooth(method = "lm",se=F)+theme_bw(base_size = 22)+scale_y_continuous(name=expression(paste("Weber's Length (nL)\n",italic("foragers"))), breaks = seq(0, 3, 0.5))+scale_x_continuous(name = "Weber's Length (mm)",breaks = seq(0.6,1,0.1))+scale_color_brewer(palette = "Dark2", name = "Raided", labels=c("No","Yes"))+ggtitle("B")+theme(plot.margin = margin(1, 1, 1, 1, "cm"))
ggsave('weblengthforagers.svg',plot=weblengthforagers,dpi=900,width =25, height = 20, units = "cm",device='svg')
## `geom_smooth()` using formula = 'y ~ x'
require(gridExtra)
## Loading required package: gridExtra
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
together<-grid.arrange(weblengthworkers, weblengthforagers, ncol=2)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
ggsave('weblengthboth.svg',plot=together,dpi=900,width =50, height = 20, units = "cm",device='svg')
This block is for worker venom volume, P1a. First, lets represent graphically.
hist(workerdata$Venom.Volume)
#A bit skewed but mostly pretty normal
venombysize = ggplot(data=workerdata,aes(x=Webers.Length.mm, y=Venom.Volume)) + geom_point() +
xlab("Weber's Length (mm)")+ylab("Venom Volume")+theme_bw(base_size = 22)+geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE)
#larger ants seem to have more venom, so it's good i'm controlling for that pattern
venombydate = ggplot(data=workerdata, aes(Date.Frozen,Venom.Volume)) +
geom_point() +
geom_smooth(method=lm, se=TRUE ) +
xlab("Date Frozen ")+
ylab("Venom Volume")+
theme_bw(base_size = 22)
We also need a model to test for overall venom per capita
#Create a column with number of ants dissected in the colony
simplified1<-na.omit(workerdata)
simplified<-simplified1
simplified<-plyr::ddply(simplified, .(Nest), transform, n.Ants.Dissected = length(Ant))
simplified1<-subset(simplified,Behavior.Code=="N")
simplified2<-subset(simplified,Behavior.Code=="F")
#IMPORTANT: THESE DO NOT WORK IF DPLYR HAS BEEN LOADED BEFORE PLYR. BIG ERROR
simplified$Webers.Length.cubed<-(simplified$Webers.Length.mm)^3
simplified<- simplified %>% group_by(Nest, n.Ants.Dissected, Final.Worker.Count,Final.Larvae.Count,Final.Pupae.Count,Final.Queen.Count.Alate,Final.Queen.Count.Dealate,Final.Male.Count, Final.Egg.Count, Day.Frozen,Pair.Group,Raided ) %>%
summarize(mean.size.cubed=mean(Webers.Length.cubed), mean_venom_volume=mean(Venom.Volume), SizeVar=var(Webers.Length.cubed),ColonyVar=var(Venom.Volume))
## `summarise()` has grouped output by 'Nest', 'n.Ants.Dissected',
## 'Final.Worker.Count', 'Final.Larvae.Count', 'Final.Pupae.Count',
## 'Final.Queen.Count.Alate', 'Final.Queen.Count.Dealate', 'Final.Male.Count',
## 'Final.Egg.Count', 'Day.Frozen', 'Pair.Group'. You can override using the
## `.groups` argument.
nrow(simplified)
## [1] 28
colonymeanvenommodel<-lme(mean_venom_volume~Raided+mean.size.cubed, data=simplified, random= ~1|Pair.Group)
summary(colonymeanvenommodel)
## Linear mixed-effects model fit by REML
## Data: simplified
## AIC BIC logLik
## -12.44567 -6.351286 11.22283
##
## Random effects:
## Formula: ~1 | Pair.Group
## (Intercept) Residual
## StdDev: 0.1622871 0.08830954
##
## Fixed effects: mean_venom_volume ~ Raided + mean.size.cubed
## Value Std.Error DF t-value p-value
## (Intercept) 0.3979070 0.3353864 13 1.186414 0.2567
## RaidedY 0.0847105 0.0335927 12 2.521692 0.0268
## mean.size.cubed 1.2019398 0.6372499 12 1.886136 0.0837
## Correlation:
## (Intr) RaiddY
## RaidedY 0.062
## mean.size.cubed -0.989 -0.113
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.3385189 -0.6473787 0.1547170 0.5135755 1.4247412
##
## Number of Observations: 28
## Number of Groups: 14
plot(colonymeanvenommodel)
qqnorm(resid(colonymeanvenommodel))
hist(resid(colonymeanvenommodel))
#A reviewer asked if we need mean size for the colony included in the model. Removing it seems to slightly decrease the p value for Raided but does not change our results qualitatively. It does increase the AIC value, though.
colonymeanvenommodelNOSIZE<-lme(mean_venom_volume~Raided, data=simplified, random= ~1|Pair.Group)
summary(colonymeanvenommodelNOSIZE)
## Linear mixed-effects model fit by REML
## Data: simplified
## AIC BIC logLik
## -10.07063 -5.038247 9.035317
##
## Random effects:
## Formula: ~1 | Pair.Group
## (Intercept) Residual
## StdDev: 0.1710766 0.0921311
##
## Fixed effects: mean_venom_volume ~ Raided
## Value Std.Error DF t-value p-value
## (Intercept) 1.0235977 0.05193082 13 19.710794 0.0000
## RaidedY 0.0918653 0.03482228 13 2.638119 0.0205
## Correlation:
## (Intr)
## RaidedY -0.335
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.24495555 -0.40213520 -0.09777063 0.59080947 1.29096471
##
## Number of Observations: 28
## Number of Groups: 14
plot(colonymeanvenommodelNOSIZE)
qqnorm(resid(colonymeanvenommodelNOSIZE))
hist(resid(colonymeanvenommodelNOSIZE))
#On average, the venom per capita is higher in raided colonies
jitter <- position_jitter(width = 0.3, height = 0)
venompercapita<-ggplot(simplified, aes(x=Raided, y=mean_venom_volume,fill=Raided))+ylab("Venom per capita (nL)")+geom_boxplot()+geom_point(position=jitter)+theme_bw(base_size = 22)+scale_fill_brewer(palette = "Set2")+theme(legend.position = "None")+scale_x_discrete(labels = c("No", "Yes"))
venompercapita
ggsave('venompercapita2.svg',plot=venompercapita,dpi=900,width =9, height = 20, units = "cm",device='svg')